RDA - Genomic offset predictions

Author

Juliette Archambeau

Published

August 26, 2023

Introduction

Most analyses conducted in this document are based on Capblancq and Forester (2021) and the associated Github repository.

Data

Genomic data

Code
# Population-based allele frequencies
# ===================================
geno <- read.csv(here("data/DryadRepo/ImputedGenomicData_AlleleFrequencies_withoutmaf.csv"),
                     row.names = 1)

# SNP sets
# ========
snp_sets <- readRDS(here("outputs/list_snp_sets.rds"))

Climatic data

Code
# Set of climatic variables
# =========================
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))


# Climatic data
# =============
# we scale the past and future climatic data at the location of the populations
# with the parameters (mean and variance) of the past climatic data.
source(here("scripts/functions/generate_scaled_clim_datasets.R"))
clim_dfs <- generate_scaled_clim_datasets(clim_var,clim_ref_adj = FALSE)

Adaptive landscape

Projecting adaptive gradient(s) across space

Adaptively enriched genetic space

Following Capblancq and Forester (2021) (and the associated Github repository), for the different sets of SNPs (one set of control SNPs and two sets of putative adaptive loci), a RDA was performed with the set of SNPs as multivariate response and the set of climatic variables as explanatory variables. For the sets of candidate SNPs, the resulting RDA spaces can be considered as ‘adaptively enriched genetic spaces’ (Capblancq and Forester 2021; Capblancq et al. 2023).

Code
runRDA <- function(snp_set, clim_var, clim_df, geno){
  
# Keep the genomic data of the selected SNPs
geno <- geno[,snp_set]
  
# RDA formula  
form_rda <- paste("geno ~ ", paste(clim_var,collapse= " + "))
  
# run RDA with only the selected SNPs
rda_outliers <- rda(as.formula(form_rda), clim_df)

}
Code
snp_sets <- lapply(snp_sets, function(x) {
  
x$rda_model <- runRDA(snp_set = x$set_snps,
                      clim_var = clim_var,
                      clim_df = clim_dfs$clim_ref,
                      geno = geno)
return(x)

}) %>% setNames(names(snp_sets))

An RDA biplot can be used to visualize the relationship between the selected SNPs and the underlying climatic variables.

Code
make_RDAbiplot <- function(x) {
  
TAB_loci <- as.data.frame(scores(x$rda_model, choices=c(1:2), display="species", scaling="none"))
TAB_var <- as.data.frame(scores(x$rda_model, choices=c(1:2), display="bp"))

eigenvalues <- summary(eigenvals(x$rda_model, model = "constrained")) %>% as.data.frame()
varexp_axis1 <- eigenvalues[2,1] *100 
varexp_axis2 <- eigenvalues[2,2] *100 


ggplot() +
  geom_hline(yintercept=0, linetype="dashed", color = gray(.80), linewidth=0.6) +
  geom_vline(xintercept=0, linetype="dashed", color = gray(.80), linewidth=0.6) +
  geom_point(data = TAB_loci, aes(x=RDA1*3, y=RDA2*3), colour = '#CD24B2', size = 2, alpha = 0.8) + #"#F9A242FF"
  geom_segment(data = TAB_var, aes(xend=RDA1, yend=RDA2, x=0, y=0), colour="black", linewidth=0.15, linetype=1, arrow=arrow(length = unit(0.02, "npc"))) +
  geom_text(data = TAB_var, aes(x=1.1*RDA1, y=1.1*RDA2, label = row.names(TAB_var)), size = 3) +
  xlab(paste0("RDA 1 (",round(varexp_axis1,1),"%)")) + 
  ylab(paste0("RDA 2 (",round(varexp_axis2,1),"%)")) +
  facet_wrap(~paste0(x$set_name, " (n=",length(x$set_snps),")")) +
  guides(color=guide_legend(title="Locus type")) +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), 
        plot.background = element_blank(), 
        panel.background = element_blank(), 
        strip.text = element_text(size=11))

}
Code
rda_biplots <- lapply(snp_sets, make_RDAbiplot)

grid_RDAbiplots <- plot_grid(rda_biplots[[1]],rda_biplots[[2]],rda_biplots[[3]], 
                         nrow=1)

ggsave(here("figs/RDA/RDAbiplots_SNPsets.pdf"), grid_RDAbiplots, width = 14, height=6)


# We save a plot for the Supplementary Information
# ================================================
plot_grid(rda_biplots[[1]],rda_biplots[[3]],nrow=1) %>% # we keep only the control and candidate SNPs
  ggsave(here("figs/RDA/RDAbiplots_SNPsets_SI.pdf"), ., width = 9, height=6)


grid_RDAbiplots

Adaptive index across the landscape

We used the scores of the climatic variables along the RDA axes to calculate a genetic-based index for each pixel of the landscape. For each RDA axis, the index is calculated as follows:

\[\sum_{i=1}^{n}a_ib_i\]

where \(i\) is one of the \(n\) climatic variables used in the RDA model, \(a\) is the score of the climatic variable \(i\) along the RDA axis and \(b\) is the standardized value for the climatic variable \(i\) at the focal pixel.

When calculating based on the set of candidate SNPs, this index can be considered as an adaptive index that provides an estimate of the adaptive genetic similarity or difference of all pixels on the landscape as a function of the values of the climatic predictors at that location.

Code
source(here("scripts/functions/project_adaptive_index.R"))

ai_proj <- lapply(snp_sets, function(snp_set){
  
project_adaptive_index(clim_var=clim_var,K=2,snp_set=snp_set)
  
})

We project the index on maps to visualize the different gradients across the maritime pine range (which can be considered as adaptive gradients when based on the sets of putatively adaptive loci).

Code
# countries borders for the map
world <- ne_countries(scale = "medium", returnclass = "sf")

ai_maps <- lapply(snp_sets, function(snp_set){
  
RDA_proj <- ai_proj[[snp_set$set_code]]

RDA_proj <- lapply(RDA_proj, function(x) rasterToPoints(x))

# The adaptive index is scaled between 0 and 1
for(i in 1:length(RDA_proj)){
  RDA_proj[[i]][,3] <- (RDA_proj[[i]][,3]-min(RDA_proj[[i]][,3]))/(max(RDA_proj[[i]][,3])-min(RDA_proj[[i]][,3]))
}

# formatting the dataframes for ggplot
TAB_RDA <- as.data.frame(do.call(rbind, RDA_proj[1:2]))
colnames(TAB_RDA)[3] <- "value"
TAB_RDA$variable <- factor(c(rep("RDA1", nrow(RDA_proj[[1]])), rep("RDA2", nrow(RDA_proj[[2]]))), levels = c("RDA1","RDA2"))

# map options
point_size=2
x_limits = c(-10, 15)
y_limits = c(31, 52)
ggtitle <- snp_set$set_name

ggplot(data = TAB_RDA) + 
  geom_sf(data = world, fill="gray98") + 
  scale_x_continuous(limits = x_limits) +
  scale_y_continuous(limits = y_limits) + 
  geom_raster(aes(x = x, y = y, fill = cut(value, breaks=seq(0, 1, length.out=11), include.lowest = T))) + 
  scale_fill_viridis_d(alpha = 0.8, direction = -1, option = "A") +
  xlab("Longitude") + ylab("Latitude") +
  guides(fill=guide_legend(title="Genetic index")) +
  facet_grid(~ variable) +
  ggtitle(ggtitle) +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), plot.background = element_blank(), panel.background = element_blank(), strip.text = element_text(size=11))

})

# We save the maps for the Github repository
# ==========================================
pdf(here("figs/RDA/AdaptiveIndex_maps.pdf"), width=10,height=6)
lapply(ai_maps, function(x) x)
dev.off()

# We save the maps for the Supplementary Information
# ==================================================

# Candidate SNPs
pdf(here("figs/RDA/AdaptiveIndexMap_CandidateSNPs_SI.pdf"), width=9,height=4.9)
ai_maps[[1]] + theme(plot.title = element_blank())
dev.off()

# Control SNPs
pdf(here("figs/RDA/AdaptiveIndexMap_ControlSNPs_SI.pdf"), width=9,height=4.9)
ai_maps[[3]] + theme(plot.title = element_blank())
dev.off()

# Show AI maps
# ============
lapply(ai_maps, function(x) x)

Genomic offset predictions

Using spatial points

Predicting the genomic offset of the studied populations under future climates.

Code
# Function to calculate the RDA genetic offset for spatial points
# ===============================================================

# clim_ref = climatic data used to fit the RDA model (the climate-of-origin of the populations)

# clim_pred = climatic data used for predictions (so either the future climate of the populations, 
            # or climate of the common gardens or the NFI plots)

# weights = if NULL, the adaptive index is not weighted by the relative importance of the RDA axes in 
# explaining the variance.

# K = number of RDA axes used to calculate the genomic offset

# snp_set = list with at least the RDA models 


pred_GO_spatial_points <- function(snp_set, 
                                   K, 
                                   clim_ref, 
                                   clim_pred, 
                                   clim_var,
                                   weights = NULL, 
                                   CG=F){
  
# weights based on the variance explained by the different axes
weights <- (snp_set$rda_model$CCA$eig/sum(snp_set$rda_model$CCA$eig))[1:K]


list_AI <- lapply(list(clim_ref,clim_pred), function(df){
  
  lapply(1:K, function(i){
  
# Below we calculate the adaptive index for the axis i
# For that, we multiply the value of the variables at the location of the population
    # by the loadings of the variables along the axis i
ai <- df %>%
  dplyr::select(any_of(clim_var)) %>% 
  apply(1, function(x) sum(x * snp_set$rda_model$CCA$biplot[,i]))

if(!is.null(weights)){ai <- ai * weights[i]}
    
  }) %>% 
    setNames(paste0("RDA", 1:length(.))) %>% 
    as.data.frame()
  
})

if(CG==F){
  
eucloffset <- unlist(lapply(1:nrow(list_AI[[1]]), function(x) dist(rbind(list_AI[[1]][x,], list_AI[[2]][x,]), method = "euclidean")))

} else if(CG==T){
  
eucloffset <- lapply(1:nrow(list_AI[[2]]), function(x) as.matrix(dist(rbind(list_AI[[2]][x,],list_AI[[1]]), method = "euclidean"))[2:(nrow(list_AI[[1]])+1),1]) %>% 
  setNames(clim_pred[,1] %>% pull()) %>% 
  as.data.frame() %>% 
  cbind(clim_ref[,1])
  
  
  }

return(eucloffset)
}
Code
snp_sets <- lapply(snp_sets, function(x){
  
x$go <- lapply(names(clim_dfs$clim_pred), function(gcm){
  
pred_GO_spatial_points(snp_set = x,
                       clim_var=clim_var,
                       K = 2, # why K=2??
                       clim_ref = clim_dfs$clim_ref,
                       clim_pred = clim_dfs$clim_pred[[gcm]],
                       weights = T)

}) %>% setNames(names(clim_dfs$clim_pred))

return(x)
})

Relationship with Euclidean distance

Code
source(here("scripts/functions/make_eucli_plot.R"))

# Calculate the Euclidean climatic distance
list_dist_env <- clim_dfs$clim_pred %>% lapply(function(clim_pred){
  
Delta = clim_dfs$clim_ref %>% dplyr::select(any_of(clim_var)) - clim_pred %>% dplyr::select(any_of(clim_var)) 
dist_env = sqrt( rowSums(Delta^2) )

})

# Main gene pools (for the figures)
gps <- readRDS(here("data/GenomicData/MainGenePoolPopulations.rds")) %>%  arrange(pop)
Code
par(mfrow=c(1,2))


lapply(snp_sets, function(x) {

lapply(names(list_dist_env), function(gcm){
  
make_eucli_plot(
  X = list_dist_env[[gcm]],
  Y = x$go[[gcm]],
  colors = gps$color_main_gp_pop,
  color_names = gps$main_gp_pop,
  ylab = "GDM genomic offset",
  legend_position="topleft",
  plot_title = paste0(x$set_name," - ", gcm))

})
}) 

Code
# We generate scatter plots for the Supplementary Information.
# ============================================================

# Axis limits
# ===========
max_go <- lapply(snp_sets[c(1,3)], function(z){
  z$go %>% unlist()
}) %>% unlist() %>% max()
max_go <- max_go + 0.01

range_eucli <- list_dist_env %>% unlist() %>% range()

# Run the function
# ================
lapply(snp_sets[c(1,3)], function(set_i) {

p <- lapply(names(list_dist_env), function(gcm){
  
make_ggscatterplot(
  x = list_dist_env[[gcm]],
  y = set_i$go[[gcm]],
  title=gcm,
  range_eucli = range_eucli,
  max_go = max_go)

})

# remove y-labels to graphs in the second column
p[[2]] <- p[[2]] + ylab("")
p[[4]] <- p[[4]] + ylab("")


# remove x-labels to graphs in the second and third rows
p[[1]] <- p[[1]] + xlab("")
p[[2]] <- p[[2]] + xlab("")
p[[3]] <- p[[3]] + xlab("")

p[[6]] <- get_legend(p[[1]])

for(i in 1:5){p[[i]] <- p[[i]]  +  theme(legend.position = "none")} 


plot_grid(plotlist=p, nrow = 3) %>% 
  ggsave(here(paste0("figs/RDA/ScatterPlotEucliDistance_",set_i$set_code,".pdf")), 
         .,
         width=7,
         height=8,
         device="pdf")

})
$cand
[1] "C:/Users/jularc/Documents/GOPredEvalPinpin/GOPredEvalPinpin/figs/RDA/ScatterPlotEucliDistance_cand.pdf"

$control
[1] "C:/Users/jularc/Documents/GOPredEvalPinpin/GOPredEvalPinpin/figs/RDA/ScatterPlotEucliDistance_control.pdf"

Maps

Code
# Function to make the genomic offset maps
source(here("scripts/functions/make_go_map.R"))

# Population coordinates
pop_coord <-  readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_noADJ.rds")))[[1]]$ref_means %>% dplyr::select(pop,longitude,latitude)

# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(x) {
lapply(names(list_dist_env), function(gcm){
x$go[[gcm]]
}) %>%  unlist()
}) %>%  unlist() %>% range()
# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0


# Generate the maps for each set of SNPs and each GCM
lapply(snp_sets, function(x) {

go_maps <- lapply(names(list_dist_env), function(gcm){
  
make_go_map(dfcoord=pop_coord,
            snp_set = x,
            gcm=gcm,
            ggtitle=gcm,
            go_limits = go_limits,
            point_size = 3)

})

legend_maps  <- get_legend(go_maps[[1]])

go_maps <- lapply(go_maps, function(y) y + theme(legend.position = "none"))

go_maps$legend_maps <- legend_maps

go_maps <-plot_grid(plotlist=go_maps)

# =====================================================
# We save the figures for the Supplementary Information
# =====================================================
if(x$set_code=="cand"){ 
  ggsave(here("figs/RDA/GOMaps_PopLocations_CandidateSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")
} else if(x$set_code=="control"){
    ggsave(here("figs/RDA/GOMaps_PopLocations_ControlSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")
}


# =========
# Add title
# =========
title <- ggdraw() + 
  draw_label(
    x$set_name,
    fontface = 'bold',
    x = 0,
    hjust = 0
  ) +
  theme(plot.margin = margin(0, 0, 0, 7))

# merge title and plots
plot_grid(
  title, go_maps,
  ncol = 1,
  rel_heights = c(0.1, 1))

  })

For each GCM, we attribute the value 1 to the top five populations with the highest genomic offset and we attribute the value 0 to the other populations. We then count the number of 1 for each population, which gives the table and map below:

Code
source(here("scripts/functions/make_high_go_pop_maps.R"))

high_go_pops <- make_high_go_pop_maps(pop_coord=pop_coord,
                                      list_go = snp_sets$cand$go,
                                      ggtitle="RDA",
                                      nb_id_pop = 5) # number of selected populations

saveRDS(high_go_pops, file = here("outputs/RDA/high_go_pops.rds"))

high_go_pops[[1]] %>% kable_mydf
pop longitude latitude GFDL-ESM4 IPSL-CM6A-LR MPI-ESM1-2-HR MRI-ESM2-0 UKESM1-0-LL sum_go
SEG -8.45 42.82 0 0 0 0 0 0
CAD -6.42 43.54 0 0 0 0 0 0
ARM -6.46 43.30 0 0 0 0 0 0
SAC -8.36 42.12 0 0 0 0 0 0
LAM -6.22 43.56 0 0 0 0 0 0
SIE -6.49 43.53 0 0 0 0 0 0
PUE -6.63 43.55 0 0 0 0 0 0
ALT -6.49 43.28 0 0 0 0 0 0
CAS -6.98 43.50 0 0 0 0 0 0
OLB -0.62 40.17 0 0 0 0 0 0
MIM -1.30 44.13 0 0 0 0 0 0
PET -1.30 44.06 0 0 0 0 0 0
PLE -2.34 47.78 0 0 0 0 0 0
HOU -1.15 45.18 0 0 0 0 0 0
STJ -2.03 46.76 0 0 0 0 0 0
SAL -3.06 41.83 0 0 0 0 0 0
VER -1.09 45.55 0 0 0 0 0 0
TAM -5.02 33.60 0 0 0 0 0 0
OLO -1.83 46.57 0 0 0 0 0 0
BAY -2.88 41.52 0 0 0 0 0 0
PIE 9.04 41.97 0 0 0 0 0 0
CAR -4.28 41.17 0 0 0 0 0 0
LEI -8.96 39.78 0 0 0 1 0 1
CEN -4.49 40.28 0 0 1 0 0 1
VAL -4.31 40.52 0 0 1 0 0 1
BON -1.66 39.99 0 0 0 0 1 1
ORI -2.35 37.53 0 0 0 0 1 1
QUA -0.36 38.97 0 0 0 1 1 2
ARN -5.12 40.19 0 0 1 1 0 2
COC -4.50 41.25 1 1 0 0 0 2
CUE -4.48 41.37 1 1 0 0 0 2
PIA 9.46 42.02 1 1 0 0 0 2
MAD -5.23 35.18 1 1 1 1 1 5
COM -3.95 36.83 1 1 1 1 1 5
Code
high_go_pops[[2]]

Comparing GO predictions

We look at the correlation across the different genomic offset predictions at the location of the populations, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code
par(mfrow=c(1,3))

lapply(names(snp_sets[[1]]$go), function(gcm){
  
lapply(snp_sets, function(x) x$go[[gcm]]) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', 
           diag = FALSE,mar=c(0,0,2,0),
           title=gcm,
           number.cex=2,tl.cex=1.5)
  
})

Using rasters

Code
# Function to calculate the RDA genetic offset using rasters
# ==========================================================

# Arguments
# =========
# snp_set = list with at least the RDA models 
# clim_var = selected climatic variables
# K = number of RDA axes used to calculate the genomic offset
# gcm = Global Climate Model of the raster with future climatic data
# clim_ref_adj = TRUE or FALSE, specify whether the point estimate climatic data used to scale the rasters should be adjusted or not for elevation
# ref_period = the reference period used to calculate the adaptive index, can be 1901-1950 or 1960-1991
# method = `loadings` or `predict` 



pred_GO_rasters <- function(snp_set, 
                                   clim_var,
                                   K, 
                                   gcm,
                                   range_buffer = shapefile(here('data/Mapping/PinpinDistriEUforgen_NFIplotsBuffer10km.shp')),
                                   method = "loadings",
                                   clim_ref_adj = FALSE,
                                   ref_period = "ref_1901_1950"){
  
# Load point estimate climatic data of the reference period
if(clim_ref_adj==TRUE){adj <- "ADJ"} else {adj <- "noADJ"}  
clim_ref_pe <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]
  
# Extract scaling parameters, i.e. mean and variance  
scale_params <- lapply(clim_var, function(x){
    
    vec_var <- clim_ref_pe$ref_means[,x] %>% pull()
    
    list(mean = mean(vec_var), sd = sd(vec_var))
    
}) %>% setNames(clim_var)
  
  
# Rasters with climatic data for the reference period
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",clim_ref_pe$range[[1]],"-",clim_ref_pe$range[[2]],"_Extent-JulietteA/"))
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
ref_rasts <- raster::stack(tif_paths)
  
# Raster with future climatic data
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
fut_rasts <- raster::stack(tif_paths)

# checking that the CRS is the same for the buffer and the rasters
if(identical(crs(range_buffer),crs(ref_rasts))==FALSE){
  stop(paste0("CRS of the range buffer is not the same as the CRS of the reference period raster."))} 
if(identical(crs(range_buffer),crs(fut_rasts))==FALSE){
  stop(paste0("CRS of the range buffer is not the same as the CRS of future climate rasters."))} 

# Mask with the range if supplied
if(!is.null(range_buffer)){
  ref_rasts <- mask(ref_rasts, range_buffer)
  fut_rasts <- mask(fut_rasts, range_buffer)}


# extract coordinates and climatic values, and mean-center the climatic data
ref_vals <- as.data.frame(rasterToPoints(ref_rasts))
for(i in clim_var){ref_vals[,i] <- (ref_vals[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}
fut_vals <- as.data.frame(rasterToPoints(fut_rasts))
for(i in clim_var){fut_vals[,i] <- (fut_vals[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}


# Predict genomic offset for each pixel based on the loadings of the climatic variables
if(method == "loadings"){
  
    ref_proj <- list()
    fut_proj <- list()
    go_proj <- list()
  
    # Projection for each RDA axis
    for(i in 1:K){
      
      # Calculate adaptive index for the reference period
      ref_rast <- ref_rasts[[1]]
      ref_rast[!is.na(ref_rast)] <- as.vector(apply(ref_vals[,clim_var], 1, function(x) sum( x * snp_set$rda_model$CCA$biplot[,i])))
      names(ref_rast) <- paste0("RDA_ref_", as.character(i))
      ref_proj[[i]] <- ref_rast
      names(ref_proj)[i] <- paste0("RDA", as.character(i))
      
      # Calculate adaptive index under future climates
      fut_rast <- fut_rasts[[1]]
      fut_rast[!is.na(fut_rast)] <- as.vector(apply(fut_vals[,clim_var], 1, function(x) sum( x * snp_set$rda_model$CCA$biplot[,i])))
      names(fut_rast) <- paste0("RDA_fut_", as.character(i))
      fut_proj[[i]] <- fut_rast
      names(fut_proj)[i] <- paste0("RDA", as.character(i))
      
      # Calculate genetic offset based on a single RDA axis
      go_proj[[i]] <- abs(ref_proj[[i]] - fut_proj[[i]])
      names(go_proj)[i] <- paste0("RDA", as.character(i))
    }
}


  
# Predict genomic offset for each pixel based on predict.RDA
  if(method == "predict"){
    
    # Prediction with the RDA model under reference-period climates and future climates
    ref_pred <- predict(snp_set$rda_model, ref_vals[,-c(1,2)], type = "lc")
    fut_pred <- predict(snp_set$rda_model, fut_vals[,-c(1,2)], type = "lc")
    
    ref_proj <- list()
    fut_proj <- list()
    go_proj <- list()
    
    for(i in 1:K){
      
      # Calculate adaptive index for the reference period
      ref_rast <- rasterFromXYZ(data.frame(ref_vals[,c(1,2)], Z = as.vector(ref_pred[,i])), crs = crs(ref_rasts))
      names(ref_rast) <- paste0("RDA_ref_", as.character(i))
      ref_proj[[i]] <- ref_rast
      names(ref_proj)[i] <- paste0("RDA", as.character(i))
      
      # Calculate adaptive index under future climates
      fut_rast <- rasterFromXYZ(data.frame(ref_vals[,c(1,2)], Z = as.vector(fut_pred[,i])), crs = crs(ref_rasts))
      names(fut_rast) <- paste0("RDA_fut_", as.character(i))
      fut_proj[[i]] <- fut_rast
      names(fut_proj)[i] <- paste0("RDA", as.character(i))
      
      # Calculate genetic offset based on a single RDA axis
      go_proj[[i]] <- abs(ref_proj[[i]] - fut_proj[[i]])
      names(go_proj)[i] <- paste0("RDA", as.character(i))
    }
  }

  # Weights based on axis eigen values
  weights <- snp_set$rda_model$CCA$eig/sum(snp_set$rda_model$CCA$eig)
  
  # Weight the current and future adaptive indices based on the eigen values of the associated axes
  ref_proj_weighted <- do.call(cbind, lapply(1:K, function(x) rasterToPoints(ref_proj[[x]])[,-c(1,2)]))
  ref_proj_weighted <- as.data.frame(do.call(cbind, lapply(1:K, function(x) ref_proj_weighted[,x]*weights[x])))
  fut_proj_weighted <- do.call(cbind, lapply(1:K, function(x) rasterToPoints(fut_proj[[x]])[,-c(1,2)]))
  fut_proj_weighted <- as.data.frame(do.call(cbind, lapply(1:K, function(x) fut_proj_weighted[,x]*weights[x])))
  
  # Predict a global genetic offset, incorporating the K first axes weighted by their eigen values
  go_global_proj <- go_proj[[1]]
  go_global_proj[!is.na(go_global_proj)] <- unlist(lapply(1:nrow(ref_proj_weighted), function(x) dist(rbind(ref_proj_weighted[x,], fut_proj_weighted[x,]), method = "euclidean")))
  names(go_global_proj) <- "global_offset"
  
  # Return projections for current and future climates for each RDA axis, prediction of genetic offset for each RDA axis and a global genetic offset 
  return(list(ref_proj = ref_proj, fut_proj = fut_proj, go_proj = go_proj, go_global_proj = go_global_proj, weights = weights[1:K]))
}
Code
# For each snp set and each GCM, we generate a list of rasters with:
  # the projection of AI for reference_period and future climates
  # genomic offset predictions for each RDA axis 
  # genomic offset predictions incorporating the K first axes weighted by their eigen values
go_rasters <- snp_sets %>% lapply(function(snp_set){
  lapply(names(clim_dfs$clim_pred), function(gcm){
  
  pred_GO_rasters(snp_set,clim_var,K=2,gcm)}) %>% 
    setNames(names(clim_dfs$clim_pred))
  
  })

go_rasters %>% saveRDS(file=here("outputs/RDA/go_rasters.rds"))
Code
go_rasters <- readRDS(file=here("outputs/RDA/go_rasters.rds"))

# Map options
# ===========
point_size=2
x_limits = c(-10, 15)
y_limits = c(31, 52)

# get maximum and minimum genomic offset values across the different GO predictions (for each snp set and each GCM)
# so that comparing the maps will be easier
go_limits <- names(snp_sets) %>% lapply(function(x){
  
lapply(names(clim_dfs$clim_pred), function(gcm){
  
rasterToPoints(go_rasters[[x]][[gcm]]$go_global_proj) %>% as_tibble() %>% pull(global_offset)

}) %>% 
    unlist()
  
}) %>% 
  unlist() %>% 
  range()

# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0

# Mapping
# =======
go_maps <- names(snp_sets) %>% lapply(function(x){
  
go_maps <- lapply(names(clim_dfs$clim_pred), function(gcm){
  
# extract genomic offset values 
go_dfs <- rasterToPoints(go_rasters[[x]][[gcm]]$go_global_proj) %>% as_tibble()


ggplot(data = go_dfs) + 
  geom_sf(data = world, fill="gray98") + 
  scale_x_continuous(limits = x_limits) +
  scale_y_continuous(limits = y_limits) + 
  geom_raster(aes(x = x, y = y, fill = global_offset), alpha = 1) + 
  scale_fill_gradient2(low="blue", mid= "yellow", high="red", 
                       midpoint=(max(go_limits[[2]])-min(go_limits[[1]]))/2,
                       limits=go_limits,
                       name = "Genomic offset") +
  xlab("Longitude") + ylab("Latitude") +
  ggtitle(gcm) +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), 
        plot.background = element_blank(), 
        panel.background = element_blank(), 
        #legend.position = c(0.85,0.15),
        legend.box.background = element_rect(colour = "gray80"),
        strip.text = element_text(size=11))

})

legend_maps  <- get_legend(go_maps[[1]])

go_maps <- lapply(go_maps, function(y) y + theme(legend.position = "none"))

go_maps$legend_maps <- legend_maps

go_maps <-plot_grid(plotlist=go_maps)


# We save the maps for the Supplementary Information
# ==================================================
if(x=="cand"){ 
  ggsave(here("figs/RDA/GOMap_CandidateSNPs_SI.pdf"), go_maps, width=11,height=7, device="pdf")
} else if(x=="control"){
    ggsave(here("figs/RDA/GOMap_ControlSNPs_SI.pdf"), go_maps, width=11,height=7, device="pdf")
}

# Add title
# =========
title <- ggdraw() + 
  draw_label(
    snp_sets[[x]]$set_name,
    fontface = 'bold',
    x = 0,
    hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

# merge title and plots
plot_grid(title, go_maps, ncol = 1,rel_heights = c(0.1, 1))
  
})

# We save the maps for the Github repository
# ==========================================
pdf(here("figs/RDA/GOmaps.pdf"), width=14,height=9.5)
lapply(go_maps, function(x) x)
dev.off()

# Show GO maps
# ============
lapply(go_maps, function(x) x)

Corr btw predictions with rasters or spatial points

We check that the genomic offset predictions we obtained with the spatial points are highly correlated (ie similar) to those obtained with the rasters.

Comment When we extract the genomic offset values from the rasters at the location of the populations, we have missing values for some populations. Indeed, some populations are not included in the buffer of the species range that I used, based on the EUFORGEN distribution and 10km around the NFI plots. We may want to modify the buffer of the species range to include those populations!

Code
# checking correlations
names(snp_sets) %>% lapply(function(x){
  
cor_go <- lapply(names(clim_dfs$clim_pred), function(gcm){
  
go_rast <- raster::extract(go_rasters[[x]][[gcm]]$go_global_proj, pop_coord[,c("longitude","latitude")])

cor(snp_sets[[x]]$go[[gcm]],go_rast,use= "complete.obs")
  
}) %>% unlist() 
  
tibble(gcm=names(clim_dfs$clim_pred),
       cor_go=cor_go)
  
}) %>% 
  setNames(names(snp_sets)) %>% 
  bind_rows(.id="snp_set") %>% 
  kable_mydf()
snp_set gcm cor_go
cand GFDL-ESM4 0.99
cand IPSL-CM6A-LR 1.00
cand MPI-ESM1-2-HR 1.00
cand MRI-ESM2-0 1.00
cand UKESM1-0-LL 1.00
cand_corrected GFDL-ESM4 0.99
cand_corrected IPSL-CM6A-LR 1.00
cand_corrected MPI-ESM1-2-HR 1.00
cand_corrected MRI-ESM2-0 1.00
cand_corrected UKESM1-0-LL 0.99
control GFDL-ESM4 0.99
control IPSL-CM6A-LR 1.00
control MPI-ESM1-2-HR 1.00
control MRI-ESM2-0 1.00
control UKESM1-0-LL 1.00

This is ok!

Validation - NFI plots

Code
# Load the climatic data of the NFI plots.
nfi_clim <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))
  
# Keep only the climatic variables of interest and scale the climatic data
nfi_dfs <- generate_scaled_clim_datasets(clim_var, clim_ref = nfi_clim$clim_ref, clim_pred = nfi_clim$clim_survey)

# calculate the genomic offset for the NFI plots
snp_sets <- lapply(snp_sets, function(x){
  
x$go_nfi <- pred_GO_spatial_points(snp_set = x,
                                   K = 2, # why K=2??
                                   clim_var=clim_var,
                                   clim_ref = nfi_dfs$clim_ref,
                                   clim_pred = nfi_dfs$clim_pred,
                                   weights = T)
return(x)
})

# checking missing data
# lapply(snp_sets, function(x) sum(is.na(x$go_nfi)))

# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(snp_set) snp_set$go_nfi) %>%  unlist() %>% range()
# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0


# map genomic offset predictions in the NFI plots 
lapply(snp_sets, function(x) {
  
p <- make_go_map(
  dfcoord= readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))[[1]] %>% dplyr::select(contains("ude")), 
  snp_set = x,
  type="NFI",
  point_size = 0.5,
  go_limits = go_limits,
  legend_position = c(0.85,0.15),
  legend_box_background = "gray80",
  y_limits = c(35, 51))

# Figure for the SI
# =================
if(x$set_code=="cand"){ 
  p_SI <- p + theme(plot.title = element_blank())
  ggsave(here("figs/RDA/NFI_GOmap_CandidateSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf")
  
  } else if(x$set_code=="control"){
    p_SI <- p + theme(plot.title = element_blank())
    ggsave(here("figs/RDA/NFI_GOmap_ControlSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf")
    }

# Show maps in the Quarto document
# ================================
p
})

We look at the correlation across the different genomic offset predictions in the NFI plots, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code
lapply(snp_sets, function(x) x$go_nfi) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', diag = FALSE,mar=c(0,0,2,0),
               number.cex=2,tl.cex=1.5)

Validation - Common gardens

Code
cg_clim <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>%  dplyr::select(cg,any_of(clim_var))
cg_coord <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,contains("ude"))
cg_names <- unique(cg_coord$cg)

# Generate scaled climatic datasets with climatic data at the location of the populations and at the location of the common gardens
cg_dfs <- generate_scaled_clim_datasets(clim_var, clim_pred = cg_clim)
  
# Predict genomic offset of each population when transplanted in the climate of the common gardens
snp_sets <- lapply(snp_sets, function(x){
  
x$go_cg <- pred_GO_spatial_points(snp_set = x,
                                  K = 2, # why K=2??
                                  clim_var = clim_var,
                                  clim_ref = cg_dfs$clim_ref,
                                  clim_pred = cg_dfs$clim_pred,
                                  weights = T,
                                  CG=T) %>% as_tibble()
return(x)
})

# Map genomic offset predictions at the locations of the populations
go_maps_cg <- lapply(cg_names, function(cg_name){

# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(snp_set) snp_set$go_cg[[cg_name]]) %>%  unlist() %>% range()
go_limits[[1]] <- 0

p <- lapply(snp_sets, function(x) make_go_map(dfcoord=pop_coord, 
                                              snp_set = x,
                                              point_size = 3,
                                              type="CG",
                                              go_limits = go_limits,
                                              cg_name=cg_name,
                                              cg_coord=cg_coord))

plot_grid(p[[1]],p[[2]],p[[3]],nrow=1)
  
}) %>% setNames(cg_names)

pdf(here("figs/RDA/GOmaps_CGs.pdf"), width=15,height=6)
lapply(go_maps_cg, function(x) x)
dev.off()

# show maps
lapply(go_maps_cg, function(x) x)

We look at the correlation across the different genomic offset predictions in the common gardens, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code
par(mfrow=c(1,3))

lapply(cg_names, function(cg_name){

lapply(snp_sets, function(x) x$go_cg) %>% 

lapply(function(set){
  set[[cg_name]]
}) %>% 
    as_tibble() %>% 
    cor() %>% 
    corrplot(method = 'number',type = 'lower', 
             diag = FALSE,mar=c(0,0,2,0),
             title=str_to_title(cg_name),
             number.cex=2,tl.cex=1.5)

})

Let’s save the genomic offset predictions for comparison with the other methods.

Code
snp_sets %>% saveRDS(file=here("outputs/RDA/go_predictions.rds"))

Session information

Code
devtools::session_info()
─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.2.2 (2022-10-31 ucrt)
 os       Windows 10 x64 (build 19044)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  English_United Kingdom.utf8
 ctype    English_United Kingdom.utf8
 tz       Europe/London
 date     2023-08-26
 pandoc   2.19.2 @ C:/Program Files/RStudio/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package           * version  date (UTC) lib source
 assertthat          0.2.1    2019-03-21 [1] CRAN (R 4.2.2)
 backports           1.4.1    2021-12-13 [1] CRAN (R 4.2.0)
 broom               1.0.2    2022-12-15 [1] CRAN (R 4.2.2)
 cachem              1.0.6    2021-08-19 [1] CRAN (R 4.2.2)
 callr               3.7.3    2022-11-02 [1] CRAN (R 4.2.2)
 cellranger          1.1.0    2016-07-27 [1] CRAN (R 4.2.2)
 class               7.3-20   2022-01-16 [2] CRAN (R 4.2.2)
 classInt            0.4-8    2022-09-29 [1] CRAN (R 4.2.2)
 cli                 3.6.0    2023-01-09 [1] CRAN (R 4.2.2)
 cluster             2.1.4    2022-08-22 [2] CRAN (R 4.2.2)
 codetools           0.2-18   2020-11-04 [2] CRAN (R 4.2.2)
 colorspace          2.0-3    2022-02-21 [1] CRAN (R 4.2.2)
 corrplot          * 0.92     2021-11-18 [1] CRAN (R 4.2.2)
 cowplot           * 1.1.1    2020-12-30 [1] CRAN (R 4.2.2)
 crayon              1.5.2    2022-09-29 [1] CRAN (R 4.2.2)
 DBI                 1.1.3    2022-06-18 [1] CRAN (R 4.2.2)
 dbplyr              2.2.1    2022-06-27 [1] CRAN (R 4.2.2)
 devtools            2.4.5    2022-10-11 [1] CRAN (R 4.2.2)
 digest              0.6.31   2022-12-11 [1] CRAN (R 4.2.2)
 dplyr             * 1.0.10   2022-09-01 [1] CRAN (R 4.2.2)
 e1071               1.7-13   2023-02-01 [1] CRAN (R 4.2.2)
 ellipsis            0.3.2    2021-04-29 [1] CRAN (R 4.2.2)
 evaluate            0.19     2022-12-13 [1] CRAN (R 4.2.2)
 fansi               1.0.3    2022-03-24 [1] CRAN (R 4.2.2)
 farver              2.1.1    2022-07-06 [1] CRAN (R 4.2.2)
 fastmap             1.1.0    2021-01-25 [1] CRAN (R 4.2.2)
 forcats           * 0.5.2    2022-08-19 [1] CRAN (R 4.2.2)
 fs                  1.6.0    2023-01-23 [1] CRAN (R 4.2.2)
 gargle              1.2.1    2022-09-08 [1] CRAN (R 4.2.2)
 generics            0.1.3    2022-07-05 [1] CRAN (R 4.2.2)
 ggplot2           * 3.4.0    2022-11-04 [1] CRAN (R 4.2.2)
 glue                1.6.2    2022-02-24 [1] CRAN (R 4.2.2)
 googledrive         2.0.0    2021-07-08 [1] CRAN (R 4.2.2)
 googlesheets4       1.0.1    2022-08-13 [1] CRAN (R 4.2.2)
 gtable              0.3.3    2023-03-21 [1] CRAN (R 4.2.3)
 haven               2.5.1    2022-08-22 [1] CRAN (R 4.2.2)
 here              * 1.0.1    2020-12-13 [1] CRAN (R 4.2.2)
 highr               0.10     2022-12-22 [1] CRAN (R 4.2.2)
 hms                 1.1.2    2022-08-19 [1] CRAN (R 4.2.2)
 htmltools           0.5.4    2022-12-07 [1] CRAN (R 4.2.2)
 htmlwidgets         1.6.0    2022-12-15 [1] CRAN (R 4.2.2)
 httpuv              1.6.8    2023-01-12 [1] CRAN (R 4.2.2)
 httr                1.4.4    2022-08-17 [1] CRAN (R 4.2.2)
 jsonlite            1.8.4    2022-12-06 [1] CRAN (R 4.2.2)
 kableExtra        * 1.3.4    2021-02-20 [1] CRAN (R 4.2.2)
 KernSmooth          2.23-20  2021-05-03 [2] CRAN (R 4.2.2)
 knitr             * 1.41     2022-11-18 [1] CRAN (R 4.2.2)
 labeling            0.4.2    2020-10-20 [1] CRAN (R 4.2.0)
 later               1.3.0    2021-08-18 [1] CRAN (R 4.2.2)
 latex2exp         * 0.9.6    2022-11-28 [1] CRAN (R 4.2.2)
 lattice           * 0.20-45  2021-09-22 [2] CRAN (R 4.2.2)
 lifecycle           1.0.3    2022-10-07 [1] CRAN (R 4.2.2)
 lubridate           1.9.0    2022-11-06 [1] CRAN (R 4.2.2)
 magrittr          * 2.0.3    2022-03-30 [1] CRAN (R 4.2.2)
 MASS                7.3-58.1 2022-08-03 [2] CRAN (R 4.2.2)
 Matrix              1.5-1    2022-09-13 [2] CRAN (R 4.2.2)
 memoise             2.0.1    2021-11-26 [1] CRAN (R 4.2.2)
 mgcv                1.8-41   2022-10-21 [2] CRAN (R 4.2.2)
 mime                0.12     2021-09-28 [1] CRAN (R 4.2.0)
 miniUI              0.1.1.1  2018-05-18 [1] CRAN (R 4.2.2)
 modelr              0.1.10   2022-11-11 [1] CRAN (R 4.2.2)
 munsell             0.5.0    2018-06-12 [1] CRAN (R 4.2.2)
 nlme                3.1-160  2022-10-10 [2] CRAN (R 4.2.2)
 permute           * 0.9-7    2022-01-27 [1] CRAN (R 4.2.2)
 pillar              1.9.0    2023-03-22 [1] CRAN (R 4.2.3)
 pkgbuild            1.4.0    2022-11-27 [1] CRAN (R 4.2.2)
 pkgconfig           2.0.3    2019-09-22 [1] CRAN (R 4.2.2)
 pkgload             1.3.2    2022-11-16 [1] CRAN (R 4.2.2)
 prettyunits         1.1.1    2020-01-24 [1] CRAN (R 4.2.2)
 processx            3.8.0    2022-10-26 [1] CRAN (R 4.2.2)
 profvis             0.3.7    2020-11-02 [1] CRAN (R 4.2.2)
 promises            1.2.0.1  2021-02-11 [1] CRAN (R 4.2.2)
 proxy               0.4-27   2022-06-09 [1] CRAN (R 4.2.2)
 ps                  1.7.2    2022-10-26 [1] CRAN (R 4.2.2)
 purrr             * 1.0.0    2022-12-20 [1] CRAN (R 4.2.2)
 R6                  2.5.1    2021-08-19 [1] CRAN (R 4.2.2)
 ragg                1.2.4    2022-10-24 [1] CRAN (R 4.2.2)
 raster            * 3.6-11   2022-11-28 [1] CRAN (R 4.2.2)
 RColorBrewer      * 1.1-3    2022-04-03 [1] CRAN (R 4.2.0)
 Rcpp                1.0.10   2023-01-22 [1] CRAN (R 4.2.2)
 readr             * 2.1.3    2022-10-01 [1] CRAN (R 4.2.2)
 readxl              1.4.1    2022-08-17 [1] CRAN (R 4.2.2)
 remotes             2.4.2    2021-11-30 [1] CRAN (R 4.2.2)
 reprex              2.0.2    2022-08-17 [1] CRAN (R 4.2.2)
 rgdal               1.6-5    2023-03-02 [1] CRAN (R 4.2.3)
 rlang               1.1.0    2023-03-14 [1] CRAN (R 4.2.3)
 rmarkdown           2.19     2022-12-15 [1] CRAN (R 4.2.2)
 rnaturalearth     * 0.3.2    2023-01-23 [1] CRAN (R 4.2.3)
 rnaturalearthdata * 0.1.0    2017-02-21 [1] CRAN (R 4.2.3)
 rprojroot           2.0.3    2022-04-02 [1] CRAN (R 4.2.2)
 rstudioapi          0.14     2022-08-22 [1] CRAN (R 4.2.2)
 rvest               1.0.3    2022-08-19 [1] CRAN (R 4.2.2)
 scales              1.2.1    2022-08-20 [1] CRAN (R 4.2.2)
 sessioninfo         1.2.2    2021-12-06 [1] CRAN (R 4.2.2)
 sf                * 1.0-9    2022-11-08 [1] CRAN (R 4.2.2)
 shiny               1.7.4    2022-12-15 [1] CRAN (R 4.2.2)
 sp                * 1.5-1    2022-11-07 [1] CRAN (R 4.2.2)
 stringi             1.7.8    2022-07-11 [1] CRAN (R 4.2.1)
 stringr           * 1.5.0    2022-12-02 [1] CRAN (R 4.2.2)
 svglite             2.1.0    2022-02-03 [1] CRAN (R 4.2.2)
 systemfonts         1.0.4    2022-02-11 [1] CRAN (R 4.2.2)
 terra               1.6-47   2022-12-02 [1] CRAN (R 4.2.2)
 textshaping         0.3.6    2021-10-13 [1] CRAN (R 4.2.2)
 tibble            * 3.1.8    2022-07-22 [1] CRAN (R 4.2.2)
 tidyr             * 1.2.1    2022-09-08 [1] CRAN (R 4.2.2)
 tidyselect          1.2.0    2022-10-10 [1] CRAN (R 4.2.2)
 tidyverse         * 1.3.2    2022-07-18 [1] CRAN (R 4.2.2)
 timechange          0.1.1    2022-11-04 [1] CRAN (R 4.2.2)
 tzdb                0.3.0    2022-03-28 [1] CRAN (R 4.2.2)
 units               0.8-1    2022-12-10 [1] CRAN (R 4.2.2)
 urlchecker          1.0.1    2021-11-30 [1] CRAN (R 4.2.2)
 usethis             2.1.6    2022-05-25 [1] CRAN (R 4.2.2)
 utf8                1.2.2    2021-07-24 [1] CRAN (R 4.2.2)
 vctrs               0.5.1    2022-11-16 [1] CRAN (R 4.2.2)
 vegan             * 2.6-4    2022-10-11 [1] CRAN (R 4.2.2)
 viridisLite         0.4.2    2023-05-02 [1] CRAN (R 4.2.3)
 webshot             0.5.4    2022-09-26 [1] CRAN (R 4.2.2)
 withr               2.5.0    2022-03-03 [1] CRAN (R 4.2.2)
 xfun                0.36     2022-12-21 [1] CRAN (R 4.2.2)
 xml2                1.3.3    2021-11-30 [1] CRAN (R 4.2.2)
 xtable              1.8-4    2019-04-21 [1] CRAN (R 4.2.2)
 yaml                2.3.6    2022-10-18 [1] CRAN (R 4.2.2)

 [1] C:/Users/jularc/AppData/Local/R/win-library/4.2
 [2] C:/Program Files/R/R-Install/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Capblancq, Thibaut, and Brenna R Forester. 2021. “Redundancy Analysis (RDA): A Swiss Army Knife for Landscape Genomics.” Methods in Ecology and Evolution.
Capblancq, Thibaut, Susanne Lachmuth, Matthew C Fitzpatrick, and Stephen R Keller. 2023. “From Common Gardens to Candidate Genes: Exploring Local Adaptation to Climate in Red Spruce.” New Phytologist 237 (5): 1590–1605.